{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helmholtz equation\n",
    "\n",
    "We present a differentiable solver for the Helmholtz equation, based on `paddle-harmonics`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import paddle\n",
    "import paddle.nn as nn\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "from math import ceil\n",
    "\n",
    "import sys\n",
    "\n",
    "from paddle_harmonics.sht import *\n",
    "from paddle_harmonics.examples import SphereSolver\n",
    "\n",
    "cmap='twilight_shifted'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = paddle.device.set_device('gpu' if paddle.device.cuda.device_count() > 0 else 'cpu')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a shallow water solver class in `shallow_water_equations.py`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize parameters:\n",
    "nlat = 128\n",
    "nlon = 2*nlat\n",
    "lmax = ceil(nlon/4)\n",
    "mmax = lmax\n",
    "\n",
    "# timestepping\n",
    "dt = 0.001\n",
    "maxiter = int(10 / dt)\n",
    "\n",
    "# initialize solver class\n",
    "solver = SphereSolver(nlat, nlon, dt, lmax=lmax, mmax=mmax, coeff=0.01).to(device)\n",
    "\n",
    "lons = solver.lons\n",
    "lats = solver.lats\n",
    "\n",
    "jj, ii = paddle.triu_indices(lmax, mmax)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "preparing the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uspec0 = solver.randspec()\n",
    "fig = plt.figure(figsize=(8, 6), dpi=72)\n",
    "solver.plot_specdata(uspec0, fig, cmap=cmap)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dudtspec = paddle.zeros([3, solver.lmax, solver.mmax], dtype=paddle.complex128)\n",
    "nnew = 0\n",
    "nnow = 1\n",
    "nold = 2\n",
    "\n",
    "uspec = uspec0.clone().to(device)\n",
    "# pde = 'allen-cahn'\n",
    "pde = 'ginzburg-landau'\n",
    "\n",
    "\n",
    "# prepare figure for animation\n",
    "fig = plt.figure(figsize=(8, 6), dpi=72)\n",
    "moviewriter = animation.writers['pillow'](fps=20)\n",
    "moviewriter.setup(fig, pde+'.gif', dpi=72)\n",
    "\n",
    "\n",
    "with paddle.no_grad():\n",
    "    for iter in range(maxiter+1):\n",
    "        t = iter*dt\n",
    "\n",
    "        if iter % 20 == 0:\n",
    "            print(f\"t={t:.2f}\")\n",
    "            plt.clf()\n",
    "            solver.plot_specdata(uspec, fig, cmap=cmap)\n",
    "            plt.draw()\n",
    "            moviewriter.grab_frame()\n",
    "\n",
    "        dudtspec[nnew] = solver.dudtspec(uspec, pde=pde)\n",
    "        \n",
    "        # forward euler, then 2nd-order adams-bashforth time steps to start.\n",
    "        if iter == 0:\n",
    "            dudtspec[nnow] = dudtspec[nnew]\n",
    "            dudtspec[nold] = dudtspec[nnew]\n",
    "        elif iter == 1:\n",
    "            dudtspec[nold] = dudtspec[nnew]\n",
    "\n",
    "        uspec = uspec + solver.dt*( (23./12.) * dudtspec[nnew] - (16./12.) * dudtspec[nnow] + (5./12.) * dudtspec[nold] )\n",
    "\n",
    "        # implicit hyperdiffusion for vort and div.\n",
    "        # uspec[1:] = swe_solver.hyperdiff * uspec[1:]\n",
    "        # switch indices, do next time step.\n",
    "        nsav1 = nnew\n",
    "        nsav2 = nnow\n",
    "        nnew = nold\n",
    "        nnow = nsav1\n",
    "        nold = nsav2\n",
    "\n",
    "\n",
    "moviewriter.finish()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
